I definitely, like the overall scheme, but one organizational thought. I'm concerned about how astropy/nddata/convolve/__init__.py works. As it stands now, if I do from astropy.nddata.convolve import convolve, I get a function, even though there's a module astropy/nddata/convolve/convolve.py. Right now this isn't necessarily a big deal because there's only one function in this module, but it still strikes me as quite confusing, and might get much worse if more functions are added later. In addition, requireing two nested layers of imports starts getting confusing when the user-level only involves a single function (and in the future, probably only a few more).

Instead, how about leaving this __init__.py empty, and moving what's in it to astropy/nddata/__init__.py, and changing the name of the package from "convolve" to "convolution" or something like that (e.g., in astropy/nddata/__init__.py, the import would be from .convolution.convolve import convolve. Then nothing is masked, and the convolution still lives in the ndddata module.

So if I understand correctly what you are saying, the user would then import it with

from astropy.nddata import convolve

This does make sense, since scipy has the same:

from scipy.ndimage import convolve

On a side note, it would be possible to add functions such as gaussian_filter() that would again mirror the scipy capabilities, but with support NaN values, and these could then be imported with:

from astropy.nddata import gaussian_filter

Anyway, I'll implement your suggestion, but I will wait for other comments before doing so.

@astrofrog - yep, that's what I had in mind, as long as the function convolve doesn't overshadow the package convolve (and that would be fixed by changing the package name to convolution). I think in general we want the most user-friendly stuff in the first-level subpackage, but always leave the deeper subpackages/modules accessible.

I've implemented all the comments! Is there anything else you can think of before I extend this to more dimensions? How many dimensions should this be extended to? I'm not sure if many people would need 4D convolution and above, but 1D, 2D, and 3D certainly seem useful.

I having a little issue with optimization. On line 71 of boundary_none.pyx, I want val * ker to be equal to zero if val is NaN and ker is zero (instead of NaN, which is the normal result). However, I can't figure out a way to do this without doubling the runtime. For example, doing:

if ker != 0.:
    top += val * ker
    bot += ker

literally doubles the runtime speed. Does anyone have any idea how to do this differently? Is there an efficient way to override the multiplication (or use a custom mult) function that will have the behavior NaN * 0 = 0?

If I try a different if statement such as:

if not isnan(val):
    top += val * ker
    bot += ker

then I don't get much of a slowdown. isnan was defined as a macro with:

cdef extern from "math.h":
    bint isnan(double x)

The example I'm using has no zeros nor nans, so in both cases the if statement always returns zero. But ker != 0. seems to take much longer than isnan(val). Any ideas on how to speed up ker != 0.?

I've found a way around my previous question about optimization which seems to work fine. I've now added some basic documentation. The next step will be to add more tests, extend to 1 and 3-dimensions, and add examples to the documentation.

I agree that 1D convolution should be the next priority, followed by 3D. I can't think of any particular reason why you'd want to implement 4D (for that matter, offhand I can't think of any common cases for 3D, but I imagine those do exist).

I'm also getting 3 test failures for your current version in OS X 10.6 (py2.7.2 32-bit): http://paste.pocoo.org/show/551273 http://paste.pocoo.org/show/551274 http://paste.pocoo.org/show/551275

@astrofrog - you may know this already, but one very useful trick for optimizing Cython code: do cython -a foo.pyx - that will generate a file foo.html that has the .pyx file, but where if you click on a line, it will show what C code that expands into. More importantly, it shows lines in yellow that are expanding into slow C code... so it shows where you might have forgotten to declare some cdef or something... Or it will reveal places where you can give compiler directives to speed things up (that's what made me realize the comments I made above).

No, this is not ready yet - I need to fix the normalization as you suggested, and implement 1-D and 3-D convolution. I was away for the last week, but should be able to make progress on this soon.

@eteq - you can now pass nested lists as input, integer values will work, and you can normalize the kernel on the fly (though default is not to).

I've implemented 1-d and 3-d convolution, and I'm now going to improve the docs a little before a final review.

This pull request is now ready for final review!

I think it might be good to organize the documentation you added slightly differently (to match other packages), but we can address that later after the content itself is merged in. (And overall, I like the doc content a lot!)

All tests pass for me, so as far as I'm concerned, you can merge when ready.

Agree about the documentation structure - I just wanted to have something there, but it needs to be integrated better into the rest of the docs.

@taldcroft, @iguananaut, @mdboom - I plan to merge this at the start of the week, in case you want to review this too.

I'm a latecomer to this thread, but I have an alternative solution to the convolution-ignoring-nans problem implemented here: http://code.google.com/p/agpy/source/browse/trunk/AG_fft_tools/convolve_nd.py

It uses FFTs (either numpy's or FFTW3's, if FFTW3 is installed) and gets around NaNs by setting them to zero, then creating a "weight" array that is 1 where the image is not nan and 0 where it is nan. The weight array is smoothed with the same kernel as the image, then divided out. In any location where the kernel does not encounter NaNs, the weight stays 1, but if there is a nearby NaN, the weight will be decreased, so the average will be over fewer pixels.

The implementation posted above works in N dimensions. It could use more & better unit tests, and it would be especially good to compare directly to the "convol" approach implemented here. I think it would be good to include both implementations, but the pull request is obviously more mature in terms of astropy compliance.

@keflavich Have you done any performance testing? I suspect this PR is faster because it's all Cython-optimized... and I'm also concerned about the memory cost of an FFT-based approach.

Is there some other advantage to the FFT approach aside from allowing arbitrary dimensionality?

I'm inclined to suggest that we merge this pull request now, and then @keflavich, if you want to re-write your code into an astropy-compliant form, you can issue a new PR later that either adds an n-dimensional option or replaces this one once it's ready (if there's a compelling reason to replace this one). Does that sound good?

I was under the impression that FFTs are generally the fastest way to
perform convolutions. They increase execution time as n log n, while
a straight convolution goes as n^2. So, the answer will depend on
image size - for small images, I would expect this pull request to be
faster, but for large images ffts should win.

The memory cost is certainly an issue. Again, I like having both
options. My version of the FFT convolution can be extra memory
intensive in favor of speed (it allows padding to the nearest
2^n+3^n+5^n dimensional shape).

While cython optimization is probably pretty fast, numpy's ffts and
fftws are really intended to be fast. My comparison of the different
ffts is here:
http://code.google.com/p/agpy/source/browse/trunk/tests/test_ffts.py

I haven't written tests for convolution yet... I'd like to see that
comparison though.

On Wed, Mar 21, 2012 at 4:14 PM, Erik Tollerud
reply@reply.github.com
wrote:

@keflavich Have you done any performance testing?  I suspect this PR is faster because it's all Cython-optimized... and I'm also concerned about the memory cost of an FFT-based approach.

Is there some other advantage to the FFT approach aside from allowing arbitrary dimensionality?

I'm inclined to suggest that we merge this pull request now, and then @keflavich, if you want to re-write your code into an astropy-compliant form, you can issue a new PR later that either adds an n-dimensional option or replaces this one once it's ready (if there's a compelling reason to replace this one). Does that sound good?

---

Reply to this email directly or view it on GitHub:
#155 (comment)

Adam

@keflavich I see your point here on the O(n log n) scaling, so I can definitely see the merits of having both available. I'm mainly concerned about the potential for confusion in having a variety of convolution options. But I think this could be easily remedied by an explicit enough naming scheme - e.g., where this pull request centers around the driver function just named convolve, there counld be a second one that has a driver function convolvefft and make sure to include a note in the docstrings indicating basically "use convolve when you are running into memory issues, and convolvefft to better scale with size."

Would you be fine with us merging this one as is, and have you submit a pull request later for the fft-based version (after you've got the doscstrings and code to astropy standards)?

Absolutely, that's essentially what I intended, I just wanted to bring
this to everyone's attention. I agree that "convolvefft" is a good
naming scheme, and that's a good way to distinguish the two.

On Wed, Mar 21, 2012 at 6:40 PM, Erik Tollerud
reply@reply.github.com
wrote:

@keflavich I see your point here on the O(n log n) scaling, so I can definitely see the merits of having both available.  I'm mainly concerned about the potential for confusion in having a variety of convolution options.  But I think this could be easily remedied by an explicit enough naming scheme - e.g., where this pull request centers around the driver function just named convolve, there counld be a second one that has a driver function convolvefft and make sure to include a note in the docstrings indicating basically "use convolve when you are running into memory issues, and convolvefft to better scale with size."

Would you be fine with us merging this one as is, and have you submit a pull request later for the fft-based version (after you've got the doscstrings and code to astropy standards)?

---

Reply to this email directly or view it on GitHub:
#155 (comment)

Adam

I also agree that we could add always add a convolvefft (or rather convolve_fft!) function after that would live side by side with the current convolve. Note also that you could try re-using some of the tests that exist for convolve, because at the end of the day, one would hope that both techniques give the same result.

I'm going to merge thie PR now since we agree that is the way to proceed.

On the convolvefft issue, I'm kind of in favor of just having a single n-dimensional convolve function in the API. There's no reason, however, that there can't be a choice of algorithms implementing it, and selectable via a simple keyword argument. Then we can have a convolution shootout over a variety of images to determine which one would make for the best default :)

But the way I see it, they're just different strategies for doing the same thing. The advantages and disadvantages of each strategy can be explained in the docstring (just as @astrofrog's implementation is already doing for the different boundary strategies).

The downside I see to this approach is that each implementation has some optional arguments that are in no way compatible with the other, which could lead to confusion. Though I tend to think that with a bit of cleaning up, and well organized documentation, this can be mitigated.

@iguananaut - perhaps you should copy this comment into #182 now that it exists as a separate PR?

Noticed while writing convolve_fft (#182) tests - convolve may change the kernel array (modify in-place). This happens if kernel.dtype.kind == 'f'. I think this can be solved by replacing the initial checking routines with "kernel = np.asarray(kernel,dtype='float')" as @mdboom suggested on #182. If you want to keep that behavior, it should be documented.

@keflavich - good catch! I'll fix this.